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r-{ Abstract. In recent years there has been considerable interest in meth- 

^~~i ods for diffeomorphic warping of images, with applications e.g. in medical 

vQ imaging and evolutionary biology. The original work generally cited is 

rvj that of the evolutionary biologist D'Arcy Wentworth Thompson, who 

demonstrated warps to deform images of one species into another. How- 
I I ever, unlike the deformations in modern methods, which are drawn from 

^^ the full set of diffeomorphism, he deliberately chose lower-dimensional 

^^ sets of transformations, such as planar conformal mappings. 

In this paper we study warps of such conformal mappings. The approach 
is to equip the infinite dimensional manifold of conformal embeddings 
with a Riemannian metric, and then use the corresponding geodesic equa- 
tion in order to obtain diffeomorphic warps. After deriving the geodesic 
equation, a numerical discretisation method is developed. Several exam- 
ples of geodesic warps are then given. We also show that the equation 
admits totally geodesic solutions corresponding to scaling and transla- 
tion, but not to affine transformations. 
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C*~) 1 Introduction 

cn 

f^ The use of diffeomorphic transformations in both image registration and shape 

Cn analysis is now common and utilised in many machine vision and image analysis 

\ I tasks. One image or shape is brought into alignment with another by deforming 

^ the image until some similarity measure, such as sum-of-squares distance between 

Ij pixels in the two images, reaches a minimum. The deformation is computed as 

%^ a geodesic curve with respect to some metric on the diffeomorphism group. 

C^ For a general treatment and an overview of the subject see the monograph by 

Younes [1] and references therein. 

The standard approach to the deformation method is to first perform an 
affine registration (principally to remove translation and rotation), and then to 
seek a geodesic warp of the image in the full set of diffeomorphisms of a fixed 
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2 Marsland, McLachlan, Modin, Perlmutter 

domain. Typically the setting is to use the right-invariant H^ metric, which 
leads to the so-called EPDiff equation (see e.g. [2]). However, in what is arguably 
the most influential demonstration of the application of warping methods - the 
evolutionary biologist D'Arcy Wentworth Thompson's seminal book 'On Growth 
and Form' [3] - Thompson warps images of one biological species into another 
using relatively simple types of transformations, so that the gross features of the 
two images match. In a recent review of his work, biologist Arthur Wallace says: 

"This theory cries out for causal explanation, which is something the great 
man eschewed. [. . .] His transformations suggest coordinated rather than piece- 
meal changes to development in the course of evolution, an issue which almost 
completely disappeared from view in the era of the 'modern synthesis ' of evolu- 
tionary theory, but which is of central importance again in the era of evo-devo. 
[. . .] All the tools are now in place to examine the mechanistic basis of trans- 
formations. Not only do we have phylogenetic systematics and evo-devo, but, 
so obvious that it is easy to forget, we have computers, and especially, in this 
context, advanced computer graphics. We owe it to the great man to put these 
three things together to investigate the mechanisms that produce the morpholog- 
ical changes that he captured so elegantly with little more than sheets of graph 
paper and, of course, a brilliant mind. " [4] 
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Table 1. Transformation groups used in some transformations in Chapter XII, 
'On the Theory of Transformations, or the Comparison of Related Forms', of [3]. 



We draw attention to two key aspects of Thompson's examples: (i) the trans- 
formations are as simple as possible to achieve what he considers a good enough 
match (see Table 1); and (ii) the classes of transformations that he considers all 
forms groups (or pseudogroups) , either finite or infinite dimensional. Mostly, he 
uses conformal transformations, a constraint he is reluctant to give up: 

"It is true that, in a mathematical sense, it is not a perfectly satisfactory 
or perfectly regular deformation, for the system is no longer isogonal; but [. . .] 
approaches to an isogonal system under certain conditions of friction or con- 
straint. " [3, p. 1064] 
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" [. . .] it will perhaps be noticed that the correspondence is not always quite ac- 
curate in small details. It could easily have been made much more accurate by giv- 
ing a slightly sinuous curvature to certain of the coordinates. But as they stand, 
the correspondence indicated is very close, and the .simplicity of the figures illus- 
trates all the better the general character of the tran.sformation." [ibid., p. 1074] 

For applications in image registration we therefore suggest to vary the group 
of transformations from which warps are drawn. If a low dimensional group gives 
a close match, then it should be preferred over a similar match from a higher- 
dimensional group. If necessary, local deformations from the full difFeomorphism 
group can be added later to account for fine details. In this paper we consider 
the case of conformal transformations. More precisely, we consider the problem 
of formulating and solving a geodesic equation on the space of conformal map- 
pings. This is a fundamental sub-task in the framework of large deformation 
diffeomorphic metric mapping (LDDMM) [5, 6, 7, 8, 9, 10, 11, 12], which is 
the standard setup for diffeomorphic image registration. Based on the geodesic 
equation derived in this paper, the full conformal image registration problem 
will be considered in future work. 

Although the composition of two conformal maps is conformal, it need not 
be invertible: we need to restrict the domain. The invertible conformal maps 
from the disk to itself do form a group, the disk-preserving Mobius group, but it 
is only 3 dimensional. We are therefore led to consider the infinite dimensional 
configuration space Con(U,M^) of conformal embeddings of a simply connected 
compact domain U C M^ into the plane. This set is not a group, but it is a 
pseudogroup. 

In [13] the authors study a geodesic equation using an L^-metric on the 
infinite dimensional manifold of conformal embeddings and a numerical method 
is developed for the initial value problem, based on the reproducing Bergman 
kernel. Using numerical examples, it is shown that the geodesic equation is ill- 
conditioned as an initial value problem, and that cusps are developed in finite 
time which leads to a break-down of the dynamics. 

In this paper we continue the study of geodesies on the manifold of confor- 
mal embeddings, but now with respect to a more general class of Sobolev type 
H^ metrics. Furthermore, we develop a new numerical algorithm for solving the 
equations. The new method is based on a discrete variational principle, and di- 
rectly solves the two point boundary value problem. Our new numerical method 
behaves well, i.e., converges fast, for all the examples we tried as long as the 
distance between the initial and final point on the manifold is not too large. 
From this observation we expect that the initial value problem with a > is 
well-posed in the H'^ Banach space topology, which would imply that the Rie- 
mann exponential is a local difFeomorphism (see [14, 15, 16]). This question will 
be investigated in detail in future work, as it is out of the scope of the current 
paper. The experimental results in this paper (Section 6) are limiting to con- 
firming that the discrete Lagrangian method can reproduce the known geodesies 
consisting of linear conformal maps and can also calculate non-linear geodesies 
with moderately large deformations. This is a first step towards exploring the 
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metric geometry of the conformal embeddings, similar to what was done for 
metrics on planar curves by Michor and Mumford [17]. 

For analysis of 2D shapes, a setting using conformal mappings is devel- 
oped in [18]. There it is shown that the space of planar shapes is isomorphic 
to the quotient space Diff(S'i)/PSL(2,M), where PSL(2,K) acts on Diff(5i) 
by right composition of its corresponding disk preserving Mobius transforma- 
tion restricted to S^. Furthermore it is shown that there is a natural metric on 
Diff(S'^)/PSL(2.IR), the Weil-Peterson metric, which has non-positive sectional 
curvature. The setting in this paper is related but different: rather than study- 
ing planar shapes we study conformal transformations between planar domains 
and we think of the manifold of conformal embeddings as a submanifold of the 
full space of planar embeddings. The equation we obtain can be seen as a gen- 
eralisation and a restriction of the EPDiff equation. First, a generalisation by 
going from the group of diffeomorphisms of a fixed domain to the manifold of 
embeddings from a planar domain into the entire plane. '^ ^ Second, a restriction 
by restricting to conformal embeddings. The approach we take is similar to (and 
much influenced by) the recent paper [19], in which a geometric framework for 
moving boundary continuum equations in physics is developed. 

2 Mathematical Setting and Choice of Metric 

The linear space of smooth maps U — >■ M^ is denoted C°°(U,M^). Recall that 
this space is a Frechet space, i.e., it has a topology defined by a countable 
set of semi-norms (see [20, Sect. 1.1] for details on the Frechet topology used). 
The full set of embeddings U — ?► K^, denoted Emb(U,M^), is a open subset of 
C°°(U,M^). In particular, this implies that Emb(U,M^) is a Frechet manifold 
(see [20, Sect. 1.4.1]). 

Since U C K^ it holds that Emb(U,M^) contains the identity mapping on U, 
which we denote by Id. The tangent space TidEmb(U,R^) at the identity is 
given by the smooth vector fields on U, which we denote by X(U). Notice that 
the vector fields need not be tangential to the boundary dU. Also notice that 
X(U) is a Frechet Lie algebra with bracket given by minus the Lie derivative on 
vector fields, i.e., if ^,77 G X(U), then ad^{ri) — ~£^ri. 

Let g = (lx®Ax + Ay®Ayhe the standard Euclidean metric on M? . Consider 
the subspace of C°°(U,E^) consisting of maps that preserve the metric up to 
multiplication with elements in the space 9^(U) of smooth real valued functions 
on U. That is, the subspace 

C,°°(U,m2) ^{cf)<E C°°(U,M2);0*g = Pg,F e J(U)}. 



3 



4 



This generalisation of EPDiff has not yet been worked out in detail in the literature. 
However, it is likely that the approach developed in [19] for free boundary flow can 
be used with only minor modifications. 

One can also look at the generalisation of EPDiff to embeddings from a Klein 
geometry perspective. Indeed, let Diffu(IR^) denote the diffeomorphisms that leaves 
the domain U invariant. Then the embeddings Emb(U,K^) can be identified with 
the space of co-sets Difr(]R^)/Difru(K^). 
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This subspace is topologically closed in C°°(U, M.'^). The set of conformal embed- 
dings Con(U,M2) = Emb{U,K2) n C;?°(U,M2) is an open subset of C^{\J,R'^) 
and a Frechet submanifold of Emb(U,M^). The tangent space TidCon(U,M^) is 
given by 

Xe(U) = UeX(U);£5g = div(0g}, 

which follows by straightforward calculations. Notice that Xc(U) is a subalgebra 
of X(U), since 

"f i'jrjg = i^c-f >)§ ~ £v£ie. = ■^c(div(77)g) - i:,,(div(^)g) 

= {£^ div(77) - £,, div(^))g + div{T])£^g - div(^)£,,g 

^ ^ ^ 



= div(i:^r/)g. 

In the forthcoming, we identify the plane M^ with the complex numbers C 
through {x,y) i-^ z ^ X + iy. Hence, the vector fields X(U) are identified with 
smooth complex valued functions on U, and Xc(U) with the space of holomorphic 
functions. 

The complex L^ inner product on X(U) is given by 



U,i))LHu)--^ / i{z)v{z)dA{z), 
Ju 

where dA = dx A dy is the canonical volume form on M^. Correspondingly, we 
also have the real L^ inner product given by 



{^,v)lhu)-= / g(e,?y)dA = Re((e,r,))i2(u). 

Also, we have the more general class of real and complex H^ inner products 
given by 

where a > and £,x,S.y respectively denotes derivatives with respect to the 
Cartesian coordinates {x,y). Notice that if ^,?7 G Xc(U), then 

((C,'7))Hi(u) = U,v))l^v) +a((C','7'»L2(u) 

where ^' and ij' denote complex derivatives. 

The class of inner products ((•, •))//! (u) on TidEmb(U,M^) = X(U) induces a 
corresponding class of Riemannian metrics on Emb(U,K^) by 

r^Emb(U,M2)xT^Emb(U,M2)9 {U,V) ^ {U o ^-\V o ^-^)hi^^^^u))- (1) 

Note that ip^^ is well-defined as a map ip{D) — > U since (p is an embedding. 
Also note that the restriction of the metric (1) to the submanifold of diffeomor- 
phisms Diff(U) C Emb(U,M^) yields to the "ordinary" H^ metric on Diff(U) 
corresponding to the EPDiff equations. 



6 Marsland, McLachlan, Modin, Perlmutter 

3 Derivation of the Geodesic Equation 

In this section wc derive the geodesic equation on Con(U,M^) for the class of 
H^{U) metrics given by (1). These equation are given by the Euler-Lagrange 
equations with respect to the quadratic Lagrangian on Con(U,M^) given by 

where (p € T<^Enib(U,M^) corresponds to the time derivative. 
As a first step, we have the foUowing result. 

Lemma 1. For any {(p,(p) € rEmb(U,M^) it holds that 

Proof. Let / be any complex valued function on U. By a change of variables 
w = (p{z) we obtain 

/ o ^-\w)dA{w) - / f{z)\^'{z)\MA{z). 



For the first term we take f{z) = ip{z)ip{z) which yields 

{{(po(p--\ipoip-^))L2^^((})) = / ip{z)ip{z)\(p'{z)\'^dA{z) = {{ip'ip,(p'(p))L2^uy 

Ju 
For the second term we take first notice that 



ly. ,.^ _(p' oip ^{w) 



if' O If ^{w) 

and then we take f(z) — ip' (z) / ip' (z) . The result now follows. D 

We are now ready to derive the Euler-Lagrange equation from the variational 
principle. Indeed, we look for at curve ip : [0, 1] — >■ Con(U,M'^) such that 

d_ r'l 

de 



I r.i'fe {t) O ipe (t) , ipe (i) O ip^ {t) )) ^1 (^(U)) dt = 

e=0 Jq Z 



for all variations (p^ : [0, 1] -> Con(U,M^) such that (pe(0) — 'p>{0), 'p>e{l) = <p(l) 
and ip>Q = if. To simplify notation we introduce 



/ d 



e=0 
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Notice that ■0(0) = ^(1) = 0. Using Lemma 1 and the fact that differentiation 
commutes with integration we obtain 







V^^"^'^' di yBVe))L^(U) + a((^', ^ ^^c))L2(u)jdi 

( ((v'V, '4>'v + '/'V»L2(u) + ol{{lp', ijj'))Lm}) ) di 

^ ^ 

' (((</'V,V'V-¥'V+ ^('pV)))l2(u) -a((<^'>'))L2(u)jdi 

(- ((^('^V),</5»)l2(u) + (((/7V,^V-'^V))l2(u) -a((.^',f/''))L2(u; 



dt, 

where in the last two equaUties we use integration by parts over the time variable 
and the fact that "0 vanishes at the endpoints. Notice that there are now no time 
derivatives on ip. Thus, by the fundamental lemma of calculus of variations we 
can remove the time integration and thereby obtain a weak equation which must 
be fulfilled at each point in time. In order to obtain a strong formulation, we 
need also to isolate ip from spatial derivatives. The standard approach of using 
integration by parts introduces a boundary integral term. In most examples 
of calculus of variations, this boundary term either vanishes (in the case of a 
space of tangential vector fields), or it can be treated separately giving rise to 
natural boundary conditions (in the case of a space where vector fields can have 
arbitrary small compact support). However, in the case of conformal mappings, 
there is always a global dependence between interior points, and points on the 
boundary (since holomorphic functions cannot have local support). Hence, we 
need an appropriate analogue of integration by parts which avoids boundary 
integrals. For this, consider the adjoint operator of complex differentiation, i.e., 
an operator dj : Xc(U) — >■ Xc(U) such that 

((e,r/))L2(U) = {{dJ^,v))LHv), Ve,r; e Xc(U). 

Notice that dJ depends on the domain U. In the case of the unit disk U = D, 
it holds that dj^{z) — dz{z^^{z)) — 2z^(z) + z^^'(z). In the general case, this 
operator is more complicated, but can still be computed under the assumption 
that a conformal embedding U — ?> D is known (see [21]). 
Using the operator dJ we can now proceed as follows 

= -((^((/3V),V''V'))l2(u) + (('pV>V-</'''0))l2(u) -a{{if',ip'))L2(^u) 

= -((I'^T',^ + wV'»)l2(u) - (((^(/3V,V'))l2(U) 

+ {{'P''f, {i^'py - V"/''))l2(u) - a((9j(,5',?A))i2(u) 
= -{{y?V> + wV + W'v' + OLdJ_(p', 0))l2(u) + ((9j((^V), V"y5))L2(u) 
= -(((l^f + adld,)(p + 0(^V + if'^) - ^97(vpV)»)l2(u)- 

Thus, this relation must hold for all holomorphic functions "0. However, the 
expression in the first slot of the inner product is not holomorphic, so it needs 
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to be orthogonally projected back to the set of holomorphic functions. Using 
Hodge theory for manifolds with boundary, one can show that the orthogonal 
complement of Xc(U) in X(U) with respect to the real inner product (•, •)L^m) 
is given by 

Xc(U)^ - {e e X(U); e(z) - a,f - idyF + dyG + id^G, F,Ge Jo(U)}, 

where 3^o(U) — {F £ J'{[));F\g\j — 0} are the smooth functions that vanish at 
the boundary. This result is obtained in [21]. Since Xc(U)^ is invariant under 
multiplication with i, it is also the orthogonal complement with respect to the 
complex L^ inner product. 

Now we finally arrive at the strong formulation of the Euler-Lagrange equa- 
tions 

^^{A{^)iP)-^dJ{^'^)^A{ip){d.,F-idyF + dyG + id,,G), 

d,^ - 0, (3) 

F\ou - G|au = 0, 

where A{(p) = |</?'P + adj d^ is the inertia operator (self adjoint with respect to 
the L^ inner product) and where the second to last equation means that ip is con- 
strained to be holomorphic. Indeed, one may think of equation (3) as a Lagrange- 
D'Alembert equation for a system with configuration space Emb(U,M^) which, 
by Lagrangian multipliers (F, G), is constrained to the submanifold Con(U, M^). 
In the special case U = D we get 

^(yi(^)(^) - ^a,(zVV) = A{v){d^F - idyF + dyG + id.,G), 

d,^ = 0, (4) 

F\d\) = G\q\j = 0. 

3.1 Weak Geodesic Equation in the Right Reduced Variable 

It is also possible to derive the geodesic equation using the right reduced variable 
^ = ifoip~^^ as is typically done for geodesic equations on diffeomorphism groups 
with invariant metric (see e.g. [22, 23, 24]). However, there is a difference between 
the setting of embeddings and that of diffeomorphism groups, since S, is defined 
on V3(U), which is not fixed in the embedding setting. Nevertheless, the "moving 
domain" S = (^(U) simply moves along the flow, i.e., points on the boundary 
follows the flow of the vector fleld f . For details of this setting in the two cases 
of unconstrained embeddings and volume preserving embeddings, see [19]. 

Let 7g be a variation of 7 as above, and let ^^ — '^^o^~^ . Using the calculus 
of Lie derivatives, direct calculations yield 

Ce = ?) + £nS, 

e=0 



d£ 

_d_ 

{£,, div(?7)^)i2(^(u)) + a{£,', div(?7)^')i2(^(u)) 



— ^ 
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where ?7 = g- -nVe ° f^^- From the variational principle 

we now obtain the weak form of the equation in terms of the variables (c^, ^) as 



^ 

a(^', 7)' + d,£^^ + £.,e + div(r?)^')L^(7(u))) di = 0. 

Passing now to the complex inner product, we use the formulas 

g(C, div(77)0 + ig(C, div(iry)0 = 2^?? 

which yields the weak equation 

1 ^ 

((e,7? + 2A,e + 2V0>L^(7(u)) 



+ a((e', 77' + 9,(r,'0 + 2£^0)LH^m) dt = 0. (5) 

Together with the equation Lp = S,oip, this is thus a weak form of the equation (3), 
but expressed in the variables {^,S,). 



4 Totally Geodesic Submanifolds 

In this section we investigate special solutions to equation (4). The approach 
for doing so is to find a finite dimensional submanifold of Con(D,M^) such that 
solutions curves starting and ending on this submanifold actually lie on the 
submanifold. 

Recall that a submanifold N C M of a Riemannian manifold (Af , g) is totally 
geodesic with respect to (M, g) if geodesies in TV (with respect to g restricted 
to A^) are also geodesies in M . For a thorough treatment of totally geodesic 
subgroups of Diff(M) with respect to various metrics, see [24]. 

Consider now the submanifold of linear conformal transformations 

Lin(D,M2) = {(^e Con(D,M2);(y9(z) =cz,ce C}. 

Proposition 1. Lin(D,C) is totally geodesic in Con(D, C) with respect to the 
H^ metric given by (1). 
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Proof. If f n- ip{t) is a path in Lin(D, C), i.e., ip{z) = cz with c e C, then 
^ — ip o if^^ is of the form ^(2) — az with a £ C Now, let t 1— > {(p, ^) fulfill the 
variational equation (5) for each variation of the form 77(2) = bz with b € C. We 
need to show that t 1— >■ ((^, ^) then fulfills the equation for any variation of the 
form ri{z) — z^ (since the monomials span the space of holomorphic functions). 
Thus, for the first term in (5) we get 

= {{if' ■^op),ip'-{f, + 477'^ - 2^'rj) o <^))i2(D) (6) 
= |c|''((az, fez*"' + 4fc6a2:'' - 2afez''))i2(n). 

where in the first line we use the conformal change of variables formula for 
integrals. Now, since the monomials are orthogonal with respect to ((•,-))d the 
expression vanish whenever fc ^ 1. For the second term, notice that ^' = a is 
constant. Now, if 77 is constant, then all the terms 77', 77'^, £nC vanish, li -q — z^ 
with k > 2, then the second term a((f', 77' + dzi^f]') + 2»f,,^')) is of the form 



a / a{b + 3ab)kz''-^dA{z)^a\c\^a{b + 3ab) / kz'''^ 

JcO Jo 

which vanishes for every a,ci E C This concludes the proof. D 

We now derive a differential equation for the totally geodesic solutions in 
Lin(D, C) in term of the variables (c, a) corresponding to f{z) = cz and ^(z) = 
az. From 99 = ^o (^ it follows that c = ac. Next, we plug the ansatz into the weak 
equation (5), and use that b vanish at the endpoints, which yields the equations 



c = ac 

(7) 
a(2c + a) — —4a c — aa 

These equations thus gives special solutions to equation (4). We obtain that 
^(c^ + ac) — 0, so we can analytically compute these special solutions. Notice 
that if a and c are initially real, then both a and c stay real, so the even smaller 
submanifold of pure scalings is also totally geodesic. 

Fig. 1 gives a visualization of total geodesic solutions where the two end 
transformations are given first by a pure scaling and second by a pure rotation. 
Notice that within the submanifold Lin(D, C), the smaller submanifold of scalings 
is totally geodesic, as is shown in the left figure. However, the submanifold of 
rotations is not, as is shown in the right figure. 

Remark 1. By using again the weak form (5) of the governing equation one 
can further show that the submanifold of translations is not totally geodesic in 
Con(D,C). Nor is the submanifold of affinc conformal transformations. 
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c = 0.2 c = e^""' 

Fig. 1. Geodesic curve from (Pq{z) = z to (pi{z) = cz for different values of c 
and a = 0. The mesh Hues show how the unit circle evolves. Notice that the 
scaling geodesic stays a scaling (left figure), whereas the rotation geodesic picks 
up some scaling during its time evolution (right figure). 



5 Numerical Discretization 

In this section we describe a method for numerical discretization of the equa- 
tions (4). The basic idea is to obtain a spatial discretization of the phase space 
variables {(p, ip) G TCon(D, M^) by truncation of the Taylor series. Thus, we use 
a Galerkin type approach for spatial discretization. For time discretization we 
take a variational approach, using the framework of discrete mechanics (cf. [25]). 
The discrete configuration space is given by 



= {(p e Con( 



');^W 



n-l 



eC}. 



Notice that Qn is an n-dimensional submanifoldof Con(D,M^). Since Con(D,M^) 
is a open subset of the vector space of all holomorphic maps on D (in the Frechet 
topology), it holds that the discrete configuration space Qn is an open subset 
of spanc{l,z, . . . ,2:"^^} ~ C". Thus, each tangent space T^Qn is identified 
with C" by taking the coefficients of the finite Taylor series. Together with the 
restricted H^ metric, Qr, 



Let U,V e r„Con( 



is a Riemannian manifold. 
[R2),andlet(afc)^^o,(&,)^^o 



respectively be their Tay- 



lor coefficients. Then it holds that 



fe=0 



which follows since 



z%z-' 



)LH 



S,j{i+l)/7T. 



12 Marsland, McLachlan, Modin, Perlmutter 

The next step is to obtain a numerical method that approximates the geodesies. 
For this we use the variational method obtained by the discrete Lagrangian on 
Qn X Qn given by 

Ld(ipk,Vk+i) ^hL[ , ) 



'Pk + l/2 
= ^{{^'k + l/2i^k+l - ipk),y^'k+l/2i^k+l - 'Pk)))L^O) + 
^(('^'fe+l - ^'ki^'k+1 - ¥'fe)>L2(D), 

where h > is the step size. The discrete action is thus given by 

JV-l 

Ad{'fia,...,fN) ^^ Ld{ipk,Vk+i)- (8) 

k=0 

Now, a method for numerical computation of geodesies originating from the 
identity and ending at a known configuration is obtained as follows. 

Algorithm 1 Given ip £ Con(D.IR^), an approximation to the geodesic curve 
from the identity element i^oiz) = z to ip is given by the following algorithm: 

1. Set ifiN = 5'nV'j where CP„ : Con(D,M^) — > Qn is projection by truncation of 
Taylor series. 

2. Set initial guess ipk = {1 — k/N)ipo + k/Nip^ for k — 1, . . . , TV — 1. 

3. Solve the minimization problem 

min Ad{(pa,...,(pN) 

ipi,...,iPN-ieQ„ 
with a numerical non-linear numerical minimization algorithm. 

Remark 2. In practical computations we use C" instead of Q„. Thus, as a last 
step one have to check that the solution obtained fulfils that t.pk S Qn, i.e., 
'^'ki^) 7^ for z G D. For short enough geodesies, this is guaranteed by the fact 
that ^'^{z) = 1. 

Remark 3. Notice that we solve the problem as a two point boundary value 
problem. Thus, we assume that the final state ip^ is known. In future work we will 
consider the more general optimal control problem, where the final configuration 
is determined by minimimizing a functional, such as sum-of-squares of the pixel 
difference between the destination and target image. 

5.1 EfHcient Evaluation of the Discrete Action 

Evaluation of the discrete action functional (8) requires computation of the Tay- 
lor coefficients of the product v'k+ii2^^k+i ~ Vk)- The "brute force" algorithm 
for doing this requires O(n^) operations. However, it is well known that FFT 
techniques can be used to accelerate such computations. Using this, we now give 
an 0(7Vnlogn) algorithm for evaluation of the discrete action. 
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Algorithm 2 Given (po, ■ ■ ■ , (pN, o-n efficient algorithm for computing the dis- 
crete action ^^((^Oi • ■ • i Vn) is given by: 

1. For each k ^ 0, . . . , N , compute the Taylor coefficients of ip'^. This requires 
0(Nn) operations. 

2. Compute 



fe=0 

This requires 0{Nn) operations. 
3. Set af; e C^" to contain the Taylor coefficients of f'k+1/2 '^^ ^^^ fi^^^ n — 1 

elements, and then zero padded. This requires 0(Nn) operations. 
4- Set bk G C^" to contain the Taylor coefficients of {ipk+i — Vk)/h as its first 

n elements, and then zero padded. This requires 0{Nn) operations. 

5. Compute 

ak^FFT{ak), bk^FFT{bk). 

This requires 0{Nnlogn) operations. 

6. Compute component-wise multiplication 

Ck = cikbk. 

This requires (3{Nn) operations. 

7. Compute 

Ck = IFFT(cfe). 

This requires (3{Nn\ogn) operations. 

8. Let (pk be the polynomial with Taylor coefficients given by Ck . Then compute 

h ^"^ 

fc=0 

This requires 0{Nn) operations. 
Finally, the discrete action is now given by Adipo^ ■ ■ ■ , Vn) — Ai -\- A2. 

6 Experimental Results 

In this section we use the numerical method developed in the previous section 
to confirm that the discrete Lagrangian method is able to calculate geodesies 
(solutions to (4)) with moderately large deformations. The purpose of these 
experiments is to demonstrate our numerical method: we leave further investi- 
gations of conformal image registration for another paper. In all the examples we 
compute geodesies starting from the identity V'o(-z) = z and ending at some poly- 
nomial (fii (of relatively low order) . We consider both simple linear and heavily 
non-linear warpings. The simulations are carried out with two different values of 
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the parameter a to illustrate how the geodesies depend on the metric. The data 
is given in Table 2. 

Fig. 2 shows the geodesies corresponding to scaling and rotation. Notice that 
the geodesies stay in the submanifold Lin(D,IR^), as predicted by Proposition 1. 
Also notice the difference between large and small a. For small a, the scaling 
coefficient behaves (almost) like ^cf = 0, which is the solution of equation (7) 

with a = 0, while the scaling coefficient behaves (almost) like g^yCi = 0, which 
is the asymptotic solution of equation (7) as a — >■ oo. 

Fig. 3 shows the geodesies corresponding to various non-linear transforma- 
tions. Although the differences in the geodesic paths for different values of a are 
small, we notice that for higher values of a the geodesic are "more regular" at 
the cost of occupying more volume. This is especially clear in Example 3, where, 
halfway thorough the geodesic, the "bump" on the right of the shape behaves 
differently for the two values of a. 

We anticipate that the method will allow us to study the metric geometry of 
the conformal embeddings as was done by Michor and Mumford [17] for metrics 
on planar curves and to determine, for example, which conformal embeddings 
are markedly closer in one metric than another, and how the geodesic paths 
differ between different metrics and between different groups, for example, by 
passing to a smaller group (e.g. the Mobius group) or to a larger one (e.g. the 
full diffeomorphism pseudogroup for embeddings). 



7 Conclusions 

Motivated by the preference of Thompson [3] for 'simple' warps in his examples 
of how images of one species can be deformed into those of a related species, we 
have derived the geodesic equation for planar conformal diffeomorphisms using 
the H^, metric. We have chosen conformal warps as they were used by Thompson, 
and are very simple. Of course, the animals that Thompson was interested in 
are actually three-dimensional, and for any number of dimensions bigger than 2 
the set of conformal warps is rather restricted. However, our intention is to start 
with conformal warps and to continue to build progressively more complex sets 
of deformations, rather than working always in the full diffeomorphism group as 
is conventional. 

The conformal warps admit the rigid transformations of rotation and trans- 
lation as special cases, and we have shown that these linear conformal trans- 
formations are totally geodesic in the conformal warps with respect to the H^ 
metric that we have considered. 

We have also provided a numerical discretization of the geodesic equation, 
and used it to demonstrate the effects of the parameters of the conformal warps. 
In future work we will use the algorithm that we have developed here to perform 
image registration based on conformal warps using the LDDMM framework and 
apply it to images such as real examples of those drawn by Thompson. 
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Annotation 


Coefficients 


Choice of a 


Type 


Example 1 (a,b) 


Co =0 


a = 0.1 in (a) 


Scaling 




ci = 0.5 


Q = 100 in (b) 




Example 2 (a,b) 


Co =0 


a = 0.1 in (a) 


Rotation 




Ci — exp(0.47ri) 


a = 100 in (b) 




Example 4 (a,b) 


Co = 0.0185475 


a = 0.1 in (a) 


Non-linear 




ci = 0.8034225 


a = 10 in (b) 






C2 = -0.13933275 








C3 = -0.23849625 








C4 = -0.18597975 








C5 = -0.0125472 








C6 = 0.18020775 








C7 = 0.27937125 







Example 5 (a,b) 



Co = 0.00674 + 0.053125i 
ci = 0.77654 + 0.103125i 
C2 =0.109424 + 0.1031251 
C3 = -0.052777 + 0.103125i 
a = -0.115049 + 0.103125i 
C5 = -0.0409141 + 0.1031251 
C6 = 0.126201 + 0.1031251 
C7 = 0.288402 + 0.1031251 



a ■ 
Q ■ 



0.1 in (a) 
: 10 in (b) 



Non-linear 



Table 2. Data used in the various examples. The polynomial for the final con- 
formal mapping is of the form fi{z) = X]fe=o '^kz'' with n = 16. Coefficients not 
listed are zero. In all examples, we used A'^ = 20 discretisation points in time. 
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Example 1 

(a) 






(b) 






Example 2 

(a) 







(b) 







Fig. 2. In Example 1, geodesies in the H^ metric connecting the identity z i— >■ z 
to z I— > O.bz are calculated using the discrete Lagrangian method. In (a), a = 
0.1 and in (b), a = 100. Both geodesies coincide with the analytic solution to 
equation (7). In Example 2, the geodesic connecting z i— >■ z to z i~> e^'^'^'z is 
calculated, again matching the analytic solution perfectly, illustrating the effect 
of a. 
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Example 3 

(a) 







(b) 







Example 4 

(a) 







(b) 







Fig. 3. Examples 3 and 4 illustrate two geodesies in the i/^ metric on conformal 
embeddings calculated using the discrete Lagrangian method. In (a), a = 0.1 and 
in (b), a = 100. In both examples, the target difFeomorphism ifi has been chosen 
to be highly non-linear (see Table 2 for the exact data used). Little difference is 
visible to the eye between the two values of a: in Example 3, a small bump on 
the right side of the boundary behaves differently. 
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